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Following the Dirac-Frenkel time-dependent variational principle, dynamics of a one-dimensional 
Holstein polaron is probed by employing the Davydov D2 Ansatz with two sets of variational pa- 
rameters, one for each constituting particle in the exciton-phonon system, and a simplified variant 
of the Davydov Di Ansatz, also known as the D Ansatz, with an additional set of phonon displace- 
ment parameters. A close examination of variational outputs from the two trial states reveals fine 
details of the polaron structure and intricacies of dynamic exciton-phonon interactions. Superradi- 
ance coherence sizes, speeds of exciton-induced phonon wave packets, linear optical absorption, and 
polaron energy compositions are also included in the study. 
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I. INTRODUCTION 

An electron in an insulating crystal induces a local lat- 
tice distortion around itself as it is excited by a photon. 
The electron together with the locally distorted lattice 
around it can be viewed as a quasiparticle which is also 
known as a polaron. Relaxation dynamics of photoex- 
cited entities such as polarons in liquids and solids has 
recently received much interest thanks to the advent of 
the ultrafast laser spectroscopy^—. It is now commonly 
accepted that dephasing and relaxation time scales in 
condensed matter are approximately picoseconds to tens 
of picoseconds. Emerging technological capabilities to 
control femtosecond pulse durations and down-to-one- 
hertz bandwidth resolutions provide novel probes on vi- 
brational dynamics and excitation relaxation. For ex- 
ample, progress in femtosecond spectroscopic techniques 
has made it possible to observe a coherent phonon wave 
packet oscillating along an adiabatic potential surface 
associated with a self-trapped exciton in a crystal with 
strong exciton-phonon interactions^. 

Ultrafast events occur on femtosecond to picosecond 
timescales, and studies of ultrafast relaxation dynamics 
are a strategic research domain from both fundamental 
and technological points of view. It is the aim of an 
ultrafast optical experiment to provide information on 
the details of temporal evolution on a femtosecond scale, 
which, in turn, offers insights into fundamental processes 
governing the dynamics. Developments in ultrafast laser 
physics and technologies now allow studies of nonequi- 
librium carrier/exciton dynamics that is previously inac- 
cessible to traditional linear optical spectroscopy. How- 
ever, theoretical studies of polaron dynamics has not re- 
ceived much-deserved attention due to inherent difficul- 
ties in obtaining reliable solutions. Previously, a time- 
dependent form of the Merrifield-type polaron wave func- 
tion with zero crystal momentum has been employed to 
yield an approximative solution to the Schrodinger equa- 
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tion that governs the ultrafast relaxation process of a one- 
dimensional molecular chair>£r— . Results show that tem- 
poral changes of the exciton coherence size and related 
energy relaxation strongly depend on the exciton transfer 
integral, the exciton-phonon coupling strength, and the 
phonon bandwidth. The applicability of the Merrificld 
wave function, however, is restricted to the narrow-band 
regime where the electronic coupling between neighbor- 
ing molecules is sufficiently weak leaving the electron- 
phonon coupling at a dominant role. In addition, in the 
presence of off-diagonal electron-phonon interactions, the 
Merrificld Ansatz is shown to fail^. 

Beyond the Merrifield Ansatz, there exist several trial 
wave functions of increasing sophistication to describe 
the polaron state in a translationally invariant man- 
ner such as the Toyozawa Ansatz^, the Global-Local 
(GL) Ansatz formulated by Zhao and coworkers in the 
early 90 11 ' 12 , and a delocalized form of the Davydov Di 
Ansatz that has been constructed very recently^. By 
using these Ansatze, we have previously investigated the 
ground state polaron energy band and the self-trapping 
phenomenon of a static Holstein polaron. Far superior 
results have been obtained compared with the Merri- 
field Ansatz^ i 12 ' 14 . Closely related to those polaron trial 
states are the Davydov Ansatze which originated from 
the theory of "Davydov soliton." Seeking to explain stor- 
age and transport of biological energy in protein struc- 
tures, Davydov proposed in 1973 that quantum units of 
peptide vibrational energy might become "self-localized" 
through interactions with lattice phonon o 15 ' 16 . Follow- 
ing his suggestion, many related studies have been car- 
ried out on the "Davydov soliton," an essentially one- 
dimensional object that maintains dynamic integrity by 
balancing the effects of nonlinearity against those of dis- 
persion. The original Davydov Ansatze include two forms 
of varying sophistication, namely, the Di— ~— and D2 
Ansatze, with the latter being a special case of the for- 
mer. 

In this work, we employ the Davydov D2 Ansatz and a 
localized form factor of the GL Ansatz, previously known 
as the D Ansatz, to study time evolution of the Holstein 
polaron following the Dirac-Frenkel time-dependent vari- 
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ational approach 24 . Time-dependent variational param- 
eters which specify the two Ansatze are obtained from 
solving a set of coupled differential equations generated 
by the Lagrangian formalism of the Dirac-Frenkel vari- 
ation. Special attention is paid to the evolution of the 
reduced single-exciton density matrix and the propaga- 
tion of exciton and phonon amplitudes from an initial 
location, where the exciton resides at t = 0, to the entire 
aggregate. 

The paper is organized as follows. In Sec. II we in- 
troduce the model Hamiltonian and discuss the nature 
of the trial states we use for dynamics studies, which is 
followed by the procedure of the time-dependent varia- 
tion. In Sec. Ill results from our dynamics calculation 
using the two trial states are displayed and discussed. 
Conclusions are drawn in Sec. IV. 



vacuum states of the exciton and the phonon field: | G) — 
0}ox|0) p h- We confine ourselves to one-exciton subspace 
for the optically excited stated since the exciton number 
is conserved in the total Hamiltonian in Eq. (TIJ. 

The spectral density 2 ^, embodying all relevant informa- 
tion of the coupled exciton-phonon system in Eq. (|4]) can 
be written as 

CUH ee — / (V m (t)V n (0)) ph e wt dt (6) 

where V n (t) is a Heisenberg representation of the exciton- 
phonon interaction at the nth site 

V n = J29^ q e iqn (b q + bl q ) (7) 
i 



II. METHODOLOGY 

We consider a one-dimensional aggregate of N 
molecules with a periodic boundary condition. There is 
only one two-level electronic system for each molecule 
coupled linearly with the phonon field. The Holstein 
Hamiltonian for the exciton-phonon system can be writ- 
ten as 2 ^— 

H = H cx + _ff cx _ p h + -ffph (1) 

with 

Hex = - J &n (Bn+1 + B n - 1 ) (2) 



H P h = J2^ q i>lb q (3) 
i 

N 

^cx-ph- 9 q "M B n(he iqn +i>le- iqn ) (4) 

q,n— 1 

Here H ex is the Hamiltonian for a single Frenkel exciton 
band in a rigid chain, and B n (B^) is the Pauli annihila- 
tion (creation) operator of an exciton at the nth site. We 
set h = 1 and assume the nearest-neighbor exciton trans- 
fer integral J mn = J<5m,n±i- -ffph is the phonon Hamil- 
tonian where b q (b q ) is the boson annihilation (creation) 
operator of a phonon with monentum q and frequency 
uj q . For simplicity, the zero-point energy is neglected, 
^cx-ph assumes the exciton is coupled linearly with the 
phonon field in a site diagonal form. 

Because there is no exciton in the ground state, the 
Hamiltonian for the ground state is represented by 

H g = |0) cx i?ph cx(0| (5) 

where |0) ex stands for the exciton vacuum. The global 
ground state is then described as a direct product of both 



and the (• • • ) p h denotes the thermal average in the free- 
phonon basis. Note that C mn (uj) represents not only 
the time correlation but also the spatial correlation for 
the exction-phonon interactions 2 ^. Here only the case of 
zero temperature is considered. Substituting Eq. ([7]) into 
Eq. one has 

C mn {u) = gy q e^ m -^5{u - u q ) (8) 

Following Tanaka 8 , the spectral density is assumed to 
have the form 

Coo(w) = ^2g q uj^S{uj -u q ) 
i 

o a 2 

= -^^ 2 -(^o) 2 (9) 
The relaxation energy is defined by 

Er EE ^ ^Idco = J2 & q EE SU (10) 

where S is the Huang-Rhys factor—, loq is the central 
energy of the phonon band, and W is the phonon energy 
band width . Here we assume a linear dispersion phonon 
band with a linear dispersion 

U}q = u + 2W(^-h (11) 

The central frequency of the phonon band luq is taken as 
the energy unit, i.e., ujq = 1. From Eqs. (f9"|)- (fTTj) . we can 
obtain g q for a given S and W. 

The Merrificld Ansatz, also known as the small-polaron 
Ansatz 5 .i£, is a simple, effective trial state to describe 
a coupled exciton-phonon system in which the exciton 
transfer integral is small compared with other system en- 
ergy scales. There is only one set of variational param- 
eters, the phonon displacements (3^(t), or alternatively, 
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its Fourier transform /3ff(i), that characterizes the Mer- 
rifield trial state I^m! K) of crystal momentum K: 

\*vl\K) = N-V 2 Y,e iKn \n) ex 

n 

x exp[- J2^n b -Jt b - ^*-„6»J]|0)ph 

n 

x exp[- 5^C8f e-^SJ - /?f V*"S ? )]|0) ph 

(12) 

In contrast, Davydov built a theory of enveloped solitons 
based on the two trial states of localized nature 

\D X {t)) =£Vn(t)££|0>ox® |A n (t)> (13) 
n 

with 

|A„(t)> ee exp{]T[A 9 „(^ _ A*„(i)S g ]}|Q) ph (14) 
i 

and 

|£> 2 (*)) = X;V'n(t)Si|0)„®|A(*)) (15) 

with 

\\{t)) ee exp{^[A 9 (t)6t - A;(i)S 9 ]}|0> ph (16) 

9 

The Di Ansatz differs from the more widely used D2 
Ansatz. In the D2 Ansatz the phonon amplitudes are 
exciton-site independent, i.e., X qn (t) = X q (t) for all n. 
The D2 Ansatz reprsents the quantum state of each 
phonon mode by a single coherent state, and hence has at 
all times a strongly classical character. In such solutions, 
the typically pulse-shaped distribution of the amplitudes 
X q describes a lattice deformation fixed in the frame of 
the lattice, to which conforms a correspondingly pulse- 
shaped distribution of exciton amplitudes ip n . 

Since the description of phonons in the D2 Ansatz is 
unsophisticated, direct correlations between the ampli- 
tudes of phonons and exciton are neglected. In many 
cases, the D2 Ansatz is inadequate to capture the com- 
plexities of the exciton-phonon system. Therefore, we 
introduce the D Ansatz, i.e., the localized form factor of 
the GL Ansatz: 

\D{t)) = ^2i> ni (t)Blexp[^2(X n2 - p n2 . ni )bU 

~(K 2 - & a - n >n a ]|0> 
n q 

-(/3* q e^ - X*)b q }\0) (17) 



The global amplitude X q can be related to the spatial 
average of X qn , and the local amplitude (3 q e~ lqn can 
be viewed as an Ansatz for the spatial variation of X qn 
around this mean value. 

The time evolution of the photoexcited state in a 
one-dimensional molecular aggregate follows the time- 
dependent Schrodinger equation. As far as the exciton- 
phonon coupling is weak, it can be treated perturba- 
tively leading to a nonlinear exciton equation with the 
use of the relaxation superoperator. In the strong cou- 
pling case, however, the perturbative method is no longer 
valid. There are several approaches to solve the time- 
dependent Schrodinger equation. In the Hilbert space, 
the time-dependent wave function |$(t)) for the Hamil- 
tonian H is parameterized by a m (t) (to — 1, M): 

|*(t)> ee \{a m (t)}). (18) 

Assume that $(£)) satisfies the time-dependent 
Schrodinger equation 

i^mt)) = H\m)- (19) 

Explicitly putting in the Hamiltonian H of Eq. (TTJ and 
writing d\<f>(t))/dt in terms of a m (t) and their time- 
derivatives a m (t) (to = 1, M), one has 

i^ t \{a m (t)}) = H\{a m (t)}) (20) 

Projecting Eq. (l20|) onto M different states \^ m ) (to = 
1, M), one obtains M equations of motion for the pa- 
rameters a m (t): 

(* m |{a m (i)},{d m (i)}) =0 (21) 

The approach we adopt in this work is the Lagrangian 
formalism of the Dirac-Frenkel time-dependent varia- 
tional method^, a powerful technique to obtain approxi- 
mate dynamics of many-body quantum systems for which 
exact solutions often elude researchers. We formulate the 
Lagrangian L as follows 

L=<*(t)|~ -#!*(*)> (22) 

From this Lagrangian, equations of motion for a m (t) and 
their time-derivatives a m (t) (to = 1, ...,M) can then be 
obtained 

d dL dL 

It can be shown that with a properly chosen set of |\P TO ), 
Eqs. (|2"l"T) and (f2"3")l can coincide. The reader is referred 
to Appendices \K[ IBl and [Cl where details of derivation to- 
gether with these M equations for D 2 , D and T>i Ansatze 
are given, respectively. Descriptions of the numerical pro- 
cedure can be found in Appendix |Dj 
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FIG. 2: The superradiance coherence size of the polaron L £ 
for J = 0.1, W = 0.1 and S = 0.5 calculated from the D 
(solid) and D2 (dashed) Ansatze. 



FIG. 1: The reduced exciton density matrix p mn calculated 
from the D Ansatz for (a) t — 7tt/ujo; (b) t = 87t/cjo; (c) 
t — 9n/u)o; (d) t = 107t/ojo. The control parameters are 
J = 0.1, W = 0.1 and S = 0.5. 

III. RESULTS AND DISCUSSIONS 

A. Reduced Exciton Density Matrix and 
Coherence Size 

The reduced single-exciton density matrix p mn (t) can 
be obtained after solving the coupled equations of varia- 
tional parameters 

p mn {t) = Tr[p(t)BlB n ] (24) 

where p(t) = is the full density matrix at 

zero temperature, and is the total polaron wave 

function at time t after the photo excitation takes place. 
Substituting the detailed form of the D2 Ansatz into the 
polaron wave function, one obtains 

p mn = (D 2 \BlB n \D 2 ) = r m (t)M) (25) 

For the D2 Ansatz, the reduced exciton density matrix 
involves only the exciton amplitude, and does not contain 
information of the phonon manifold. In contrast, the 
more sophisticated trial state, the D Ansatz, includes 
the Debye- Waller factor S n m in the expression of the 
reduced exciton density matrix: 

p mn = (D\BlB n \D) = \* m (t)\ n (t)S m ,„ (26) 

where S n ^ m can be written as 

q 

expIiV- 1 J2(P> tqn ~ K)^' lqm ~ A «)] 

exphl^iV-^l^e^-A,! 2 ] (27) 
q 



Using the Merrifield Ansatz^, when J = 0, solutions to 
coupled equations of variational parameters can be ob- 
tained analytically for zero phonon bandwidth (W = 0). 
It is found that L p = 1/N for all t = 2ttti/ll1o (n — 
1,2,3, ...). Periodic changes of the polaron structure as 
a function of time have also been revealed when J and 
W are both small. At uj^t = 7tt, 8tt, 9tt and 10-7T, the re- 
duced exciton density matrix p mn calculated from the D 
Ansatz for the case of J = 0.1, W = 0.1 and S = 0.5 
is displayed in Fig. [T] At t = 0, p mn is confined to the 
(0, 0) point, i.e., the exciton is localized at one site of the 
one-dimensional system. Then, as the excited state re- 
laxes, the non-zero elements of p mn spread out gradually. 
Oscillatory behavior is visible during the evolution of the 
density matrix as both the J and W are small, and as a 
result, dephasing due to phonon dispersion and exciton 
transfer is at a minimal level. When t equals to 2nir/oJo, 
the matrix is rather delocalized and have a symmetric 
distribution, and at t = (2n + l)7r/wo, matrix elements 
aggregate along the diagonal line. 

When J is increased, the oscillatory behavior of the 
density matrix gradually fades as the exciton density de- 
localization happens more swiftly. From the reduced ex- 
citon density matrix, one can obtain an exciton coherence 
size LfT^— characterizing the size of a domain within 
which the chromophores emit coherently. All informa- 
tion relevant to the excitonic superradiance is contained 
within the reduced exciton density matrix p m n{t), the 
time dependence of which provides a clear view of the 
relaxation process. Previously, for translationally invari- 
ant states, a definition of the characteristic coherence size 
L p in terms of the density matrix was introduced by ap- 
plying the inverse participation ratio concept commonly 
used in the theory of quantum localization: 

L p ^[L a Y.\p mn \ 2 \- 1 [{Y.\p^ ( 28 ) 

mn run 

This definition of the coherence size is borrowed to quan- 
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tify the excitonic delocalization despite that our trial 
states are both of the localized form. Values of L p cal- 
culated from the D2 and D trial states for the case of 
J = 0.1, W = 0.1 and S = 0.5 are shown in Fig. [2] 
as a function of time. At t = 0, the coherence size is 
1/32 because the exciton is localized to one site com- 
pletely. Then coherence size will increase in an oscillatory 
fashion. When t equals 2nir/uJo, the coherence size will 
reach a local maximum and then decrease until t equals 
(2n + 1)tt/u)q. This is consistent with the time evolution 
of the reduced exciton density matrix. 



B. Exciton Amplitudes and Phonon Displacements 
of the D2 Ansatz 

At variance with the single-parameter Merrificld 
Ansatz, the localized D2 trial state employs two sets 
of variational parameters for independent descriptions 
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FIG. 3: Time evolution of variation parameters of the D2 
Ansatz for J = 0.5, W = 0.8 and S = 0.5. (a) A bird's- 
eye view of the exciton probability |'i/i n (t)| 2 in real space; (b) 
a contour plot of |V>n(£)| 2 ; (c) a contour plot of the phonon 
displacement |A n (t)| in real space. 



of the two types of constituting particles in the coupled 
exciton-phonon system. By studying the outputs of time- 
dependent variational parameters, namely, the exciton 
amplitude ip n {t) and phonon displacement A„(i), one is 
able to probe the polaron dynamics which follows photo- 
excitation at a single site of a one-dimensional molecular 
chain. At t = 0, it is assumed that the exciton is gen- 
erated on one site and there are no initial phonon dis- 
placements on the chain. Thanks to the exciton transfer 
integral and the exciton-phonon coupling in the Holstein 
Hamiltonian, one witnesses the spreading of the exci- 
ton amplitude and the growth of phonon deformations 
around the initial exciton site. Here we use three exam- 
ples to illustrate the time evolution of the exciton am- 
plitude and the phonon displacement of the D2 Ansatz 
in the site space. As shown in Figs. |31 2] and [3] for three 
sets of (J,W), (0.5,0.8), (0.5,0.1) and (1.0,0.8) respec- 
tively, the exciton amplitude propagates from the initial 
site n — to the entire chain with a speed that is propor- 




FIG. 4: Time evolution of variation parameters of the D2 
Ansatz for J = 0.5, W = 0.1 and S = 0.5. (a) A bird's- 
eye view of the exciton probability |V'n(^)| 2 m rea l space; (b) 
a contour plot of |^i„(f)| 2 ; (c) a contour plot of the phonon 
displacement |A n (i)| in real space. 
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tional to the transfer integral J. The Huang-Rhys factor 
S is set at 0.5. In Figs. [3J and U (J = 0.5), the hrst trav- 
eling wave of the exciton amplitude reaches the the op- 
posite end of the ring (n = 16) at t « 2.6(2tt/ujo), while 
in Fig. [5] ( J = 1.0), the exciton amplitude reaches the 
n = 16 site at t w 1.3(27r/wo)- This is because the prop- 
agation speed of the exciton, or equivalently, dE cyi {k) / dk 
with E cx (k) the bare exciton band, is proportional to the 
exciton transfer integral J according to Eq. ([2]). The ef- 
fect of phonon dispersion on the exciton is reflected in the 
differing degree of phonon-induced dissipation in Figs.[3jD 
and HJd. Exciton wave packets in Fig. 2}d survive much 
longer than those in Fig. [3}j thanks to a much smaller 
phonon band width. 

At t = there are no phonon deformations anywhere 
on the one-dimensional ring. Due to photo-excitation at 
t = 0, a high concentration of the exciton density at a 
single location, appears in the vicinity of n = for a 
short time duration near t — [the red-colored spot near 
the point (0, 0) in Figs. [3b, Hp and [5b], triggering a pair 
of localized phonon wave packets which travel at group 
velocities ±v q with 

2W 

v q =V q u q = . (29) 

As shown in Figs. [3J;, 0J; and [SJ;, the angle between the 
two trajectories of phonon wave packets departing from 
point (0, 0) but with opposite velocities in Figs.[3p and [5]: 
(W = 0.8) is about 8 times of that in Fig. Hfc (W = 0.1). 

For a larger value of transfer integral, e.g., J = 1.0 
in Fig. [5J the left-moving and right-moving wave packets 
of the exciton depart the site of creation and makes a 
quick rendezvous at the opposite site of the ring (n = 16), 
where the recombined exciton density remains sufficiently 
high to trigger another pair of localized phonon wave 
packets as is clearly demonstrated in Figs. [SJd and [5];. 
The exciton will travel at a much reduced speed after 
the rendezvous at n — 16, and further stimulations of 
phonon wave packets are all but impossible for J = 1.0. 

C. Vibrational Amplitude of Exciton and Phonon 
using D Ansatz 

In Figs. [BIE] and [5J time evolution of the variational pa- 
rameters \ip n (t)\ 2 , |A„(t)| and |A„(i)| is displayed for the 
D Ansatz. When J and W are both small, the phonons 
are mostly localized in the vicinity of the site where the 
exciton is initially created. The first set of phonon dis- 
placements, |A n (£)|, will extend slightly to the neighbor- 
ing sites as the time increases because of a nonzero J 
and a finite phonon dispersion (Fig. [7^). But as shown 
in Fig. [5^, there is almost no spreading of j3 n (t), sug- 
gesting that this portion of phonon displacements closely 
follows the exciton and appears almost entirely on the 
site of exciton creation. Fig. [5Ji also reveals oscillatory 
behavior of \/3 n (t)\: when t equals (2n + 1)tt/uj , \f3 n (t)\ 
reaches its maximum value; and when t equals 2nir/ujQ, 




FIG. 5: Time evolution of variation parameters of the D2 
Ansatz for J = 1.0, W = 0.8 and S = 0.5. (a) A bird's- 
eye view of the exciton probability |V'n(^)| 2 m rea l space; (b) 
a contour plot of |^n(£)| 2 ; (c) a contour plot of the phonon 
displacement |A re (i)| in real space. 



the amplitude is at its minimum. Because the transfer 
integral J is nearly zero, the exciton stays very close to 
the site of creation. From Fig. [6^,, it is obvious that the 
exciton stays localized on one site with hardly any prop- 
agation for a small value of transfer integral J = 0.1. In 
this case, it seems that oscillations of phonon amplitudes 
also influence the movement of the exciton. The propa- 
gation of the exciton from the n = site to its nearest 
neighbors occurs only at t = 2n7r/wo when the phonon 
amplitude reaches a minimum value. Otherwise, because 
of the coupling with the phonons, the exciton will be 
trapped and unable to transfer to other sites. 

As J is increased, the exciton will have the ability to 
move away from the site of creation. As shown in Figs.EJa 
and [5]i, the larger J is, the faster the exciton transfers. 
Although J has no direct influence over the phonon dis- 
placements, a large J induces more excitonic coherence 
between adjacent sites, and in turn, causes more phonon 
deformations on those sites (cf. Fig. [7]) . The propagation 
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FIG. 6: Time evolution of the exciton probability IVVjWI 2 in 
real space using D Ansatz. The unit of time is 1/loq- ( a ) J = 
0.1, W = 0.1, S = 0.5; (b) J = 0.5, W = 0.8, S = 0.5; (c) 
J = 0.5, W = 0.1, S = 0.5; (d) J = 1.0, W = 0.8, S = 0.5; 
(e) J = 1.0, W = 0.8, S = 2.0. 



FIG. 7: Time evolution of the phonon displacements |A„(t)| 
of the D Ansatz in real space. The unit of time is 1/ojo- (a) 
J = 0.1, W = 0.1, S = 0.5; (b) J = 0.5, W = 0.8, S = 0.5; 
(c) J = 0.5, W = 0.1, S = 0.5; (d) J = 1.0, W = 0.8, S = 
0.5; (e) J = 1.0, W = 0.8, S = 2.0. 



of the localized phonon wave packets related to j3 n (t) be- 
comes faster J is increased from 0.5 to 1.0 as shown in 
Fig. E 

As for the D Ansatz, an increase in the width of 
the phonon dispersion W will substantially deepen the 
phonon dissipative effect on the exciton amplitude ip n (t) , 
and as a result the exciton coherence quickly vanishes 
from the case of W = 0.8. (cf. Fig. ^jp and |U). W 
is also one main factor which determines the velocity of 
the phonon wave packets for both types of phonon dis- 



placements X n (t) and P n (t). Fig. [7J; shows that when 
W is small, most of X n (t) will be localized on the n = 
site. As the phonon band width increases, X n (t) will form 
two traveling wave packets with opposite directions and 
propagate to the entire chain as shown in Fig. [7}d. 

The dynamics of /3 n (t) is more complex. As shown in 
Figs.[7p and[7j2, the phonon displacements described by 
parameter X n (t) propagate from the n = site to the 
rest of the aggregate. But the same is not true for (3 n (t). 
Most of j3 n (t) will stay localized on the site of exciton 
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FIG. 9: Speeds of localized phonon wave packets as a function 
of the half width of the phonon dispersion W for the two sets 
of phonon displacements X n {t) (in black) and /3 n (t) (in red). 
J = 0.5, 5 = 0.5. 



early, and the slopes are found to be the same, as shown 
in Fig. M 

Lastly, we examine the effect of the Huang- Rhys factor 
S which is a dimensionless parameter representing the 
exciton-phonon coupling strength. S can be obtained 
directly from absorption and fluorescence spectra, as it 
controls the vibrational progression that accompanies an 
exciton transition. Comparing Fig. [6jl with Fig. [6b, the 
exciton is found to be much less mobile as S is increased 
from 0.5 to 2.0, undergoing the so-called self-trapping 
transition. The fact that Fig. |6ji displays a much faster- 
moving polaron than Fig. [(J is also reflected in Fig.[5bl, in 
which the localized phonon wave packet in /3 n (t) is found 
to be much longer-lasting than that in Fig. [5b, and its 
speed greater. In addition, Fig. [7ji shows a few scattered 
localized phonon wave packets of X n (t), while Fig. [7fe 
contains only a single localized phonon wave packet of 

An(t). 



FIG. 8: Time evolution of the phonon vibrational amplitude 
|/3n(t)| in real space using D Ansatz. The unit of time is 1/wo- 
(a) J = 0.1, W = 0.1, S = 0.5; (b) J = 0.5, W = 0.8, S = 
0.5; (c) J = 0.5, W = 0.1, S = 0.5; (d) J = 1.0, W = 
0.8, S = 0.5; (e) J = 1.0, W = 0.8, S = 2.0. 



creation, but we can also find a small portion of /3 n (t) 
propagating to the entire molecular chain as we increase 
W (cf. Fig. EJd and [St). From the definition of the D 
Ansatz, Eq. (fT7|) . /3 n (t) represents the portion of phonon 
displacements that ride along the exciton, while X n (t) 
labels the phonon displacements that are not explicitly 
linked to the exciton. We therefore expect that the prop- 
agation speeds of the localized phonon wave packets from 
the two types of phonon displacements X n (t) and /3 n (t) 
are different. But both of them depend on the W lin- 



D. Absorption Spectra and System Energies 

Details on how to calculate linear absorption spectra 
of a one-dimensional exciton-phonon system from the D2 
and D Ansatze can be found in Appendix [E] Fig. [T0k 
and [TUb show two examples of the spectra get by D2 
Ansatz and D Ansatz, respectively, and the results are 
compared with those from the Merrifield Ansatz. As 
shown in Fig. [T0k . when J is small (e.g., J — 0.1) and 
S is large (e.g., S = 6.0), the spectra obtained by the 
D2 and Merrifield Ansatze are almost the same. The 
phonon sidebands in Fig. [TUb is labeled by n — 0, 1, 2, ... 
from left to right. The left most sideband, n — 0, corre- 
sponding to the zero-phonon line, is located at u> = —Sujq 
[Eq. pTJ)) ]. For S = 6.0, the zero-phonon line is located at 
u = — 6wq as shown in Fig. [TUb. According to the Huang- 
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I K 1 . 10: Linear absorption spectra of a 32-site one- 
dimensional ring of a coupled exciton-phonon system, (a) 
Using the D 2 and Merrifield Ansatze for J = 0.1, W = 0.1 
and S = 6.0. The decay factor used in the Fourier transfor- 
mation is 0.01. (b) Using the D and Merrifield Ansatze for 
J = 0.4, W = 0.8 and S = 1.0. The decay factor used in the 
Fourier transformation is 0.05. 



Rhys theory 30 , the phonon sidebands at zero temperature 
follow a Poisson distribution: 



S r ' 



-Fabs(w) = exp(-S) V* — -S(oj + E r 

— ^ n 



nw j 



(30) 



n=0 



From Eq. ([3"U||. the tallest of phonon sidebands should 
be the n = S — 1 peak when S » 1 and the overall 
distribution approaches a Gaussian. Thus, for S = 6.0, 
the tallest peak is at n = 5 as shown in Fig. [TTJk . 

Fig. [TUb displays the absorption spectrum calculated 
from the D and Merrifield Ansatze for J = 0.4, W = 0.8 
and S = 1.0. There are fine features with tiny peaks on 
the one-phonon manifold on the right side of the zero- 
phonon line. These small peaks are attributed to the 
32 phonon modes with momenta q = 2?:n q /N (N = 
32; n q = — 15, — 14, 14, 15, 16). As in this paper we 
assume a linear phonon dispersion as defined in Eq. (jTTJ) , 
there are altogether 17 values of uj q for N = 32. Each 
value of Lu q should correspond to a small peak on the 
spectrum if the coupling g q [Eq. (|7J|] is not too small, 
and from Eq. (fTTj). one can obtain the locations of the 
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FIG. 11: System energies calculated from the D2 and D 
Ansatz. The results are compared with those of the Mer- 
rifield Ansatz. (a) J = 0.1, W = 0.1, S = 4.0; (b) 
J = 0.1, W = 0.8, S = 4.0; (c) J = 1.0, W = 0.1, S = 0.5; 
(d) J = 1.0, W = 0.8, S = 0.5. 



small peaks u! sp (n q ) (n q = 0, 1, N/2): 
u sp (n q ) = a; ZPL + u> - W + 



4W 
7 

N 



(31) 



where ljzpl is the location of the zero-phonon line. For 
W = 0.8w and N = 32, u sp (n q ) = cj zpl + (0.2 + 0.1 * 
n q )uJo as shown in the inset of Fig. [TOb . 

Lastly, four cases of system energies calculated from 
the D2 and D Ansatze are shown in Fig. [TT] and com- 
pared with those from the Merrifield Ansatz . Since the 
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Hamiltonian defined in Eq. ([T]) is time-independent, the 
total energy of the system is expected to be a constant 
during the time evolution. At t = 0, an exciton is as- 
sumed localized at one site and there are no phonons on 
the entire ring. As shown in Fig. 1111 the phonon energy 
i?ph rises from zero and oscillates as time goes on. In 
the meantime, the interaction energy between the exci- 
ton and phonons, _E C x-ph, oscillates with almost the same 
amplitude but an opposite sign. And the total energy 
Etot of the system stays as a constant at all times. As 
shown in Figs.lTTk andlTTb. the D2 Ansatz and Merrifield 
Ansatz agree with each other when J is small and S is 
large. However, when J is large (e.g., J = 1.0) and S is 
small (e.g., S = 0.5), these two Ansatze no longer agree. 
For this case, the Merrifield Ansatz, which is translation- 
ally invariant (a Bloch wave function) but with a built-in 
small-polaron correlation, is not suitable to describe the 
dynamics. As shown in Figs. HTb andlTTH. the total sys- 
tem energy E to t vanishes during the time evolution of the 
D Ansatz. While in the Merrifield Ansatz, E to t remains 
at a constant value —2 that is equal to the initial value 
of E ex . This is because in the Merrifield Ansatz, the ex- 
citon amplitude is assumed to distribute uniformly over 
all the sites of the system, thus according to Eq. ([2]), E ex 
has a negative value when J is not negligible. However, 
for an initial state in which the exciton is localized at one 
site, according to Eq. @, E cx should be zero at t = as 
shown in Figs.fTTb andfTTH for the D Ansatz. This shows 
that the D Ansatz is a more flexible trial state than the 
Merrifield Ansatz. 

Lastly, we give a brief discussion on the validity of 
the Davydov Ansatze. A Davydov Ansatz is an approx- 
imative solution to the Schrodinger equation with the 
Holstein Hamiltonian. For a trial wave function \D(t)) 
that does not strictly obey the Schrodinger equation, the 
deviation vector \8(t)) can be written as 

\S(t))=ih^\D(t))-H\D(t)) (32) 

Here H is the Holstein Hamiltonian Eq. (1). For the 
Davydov Di Ansatz, the explicit form of \8(t))H was 
given by Skrinjar et al. in 1988. It was also proven that 
\8(t)) is orthogonal to \Di(t)). However, such orthogonal- 
ity relations are insufficient to conclude that the devia- 
tion vector \S(t)) is negligible, and the trial state is a good 
approximation to the true solution of Schrodinger equa- 
tion. To have a quantitative measure of the Schrodinger- 
equation deviation, one needs to calculate the amplitude 
of the deviation vector \S(t)), which is defined as A(t): 



aw = vfii 



(33) 



An explicit expression for A(i) as the D2 Ansatz is sub- 
stituted into the Schrodinger equation is derived in Ap- 
pendix |FJ 

Note that the dimension of A(i) is that of the energy. 
Therefore, one can gauge whether the trial state is a good 
approximative solution by comparing A(i) with the sys- 
tem energies. Two examples of such comparisons for the 
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FIG. 12: Amplitude of the deviation vector \S(t)) calculated 
from the D2 Ansatz. The results are compared with the 
system energies, (a) J = 0.1, W = 0.1, S = 4.0; (b) 
J = 0.1, W = 0.8, 5 = 4.0. 



D2 Ansatz arc shown in Figs. [r2b and IT2"b . The control 
parameters in Figs.fT2k and [12b are the same as those in 
Figs. ITTa. andfTTb. respectively. For both two cases, the 
main system energies are E p h and -E cx _ p h, and A(t) is 
found to be negligible to either E p h or i^x^ph. We con- 
clude that the D2 Ansatz yields quantitatively accurate 
solutions to the Schrodinger equation for these two cases. 



IV. CONCLUSION 

In this paper we simulate polaronic dynamics in a one- 
dimensional molecular chain following the Dirac-Frenkel 
time-dependent variational approach. After the optical 
excitation, the coupled exciton-phonon system will un- 
dergo a relaxation process from the initial photo-induced 
nonequilibrium. Based on the Holstein Hamiltonian, we 
examined time evolution of the exciton amplitude, the 
reduced exciton density matrix, the exciton coherence 
size, linear optical absorption, and induced phonon dis- 
placements for two types of variational wave functions, 
namely, the Davydov D2 Ansatz, and a simplified variant 
of Davydov Di Ansatz, also known as the D Ansatz. It 
is shown that following the equations of motion derived 
from the time-dependent variation, the exciton ampli- 
tude will transfer from the site of creation to neighboring 
sites, and as J increases, the velocity of exciton propaga- 
tion will increase as well. Exciton-induced phonon defor- 
mations can be found to form at the locations where the 
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exciton has considerable densities, and then propagate 
to the entire aggregate if there is sufficient phonon dis- 
persion enabling mobility. The half-width of the phonon 
dispersion, W , determines the speeds of phonon wave 
packets because of a linear phonon dispersion relation 
given in the model. For a large value of the transfer 
integral, the left-moving and right-moving wave packets 
of the exciton departing from the site of creation may 
make a quick rendezvous at the opposite end of the ring, 
and therefore, trigger a second pair of localized phonon 
wave packets. Linear absorption spectra can be derived 
from the exciton amplitude and the phonon displace- 
ments computed from the D2 and D Ansatze. Various 
system energies are also calculated and analyzed for the 
one-dimensional coupled exciton-phonon molecular ring 
after photo-excitation. Overall it is found that D is a 
more flexible trial state, while the D2 Ansatz is rather effi- 
cient for computation, and its extension to higher spatial 
dimensions is a feasible generalization of our approach. 

Our trial states in this paper are both localized wave 
functions, and we intend to work out detailed dynam- 
ics of their translationally invariant counterparts, the 
Toyozawa Ansatz and the GL Ansatz. Our approaches 
here can be also readily extended to include other forms 
of exciton-phonon interactions, such as asymmetric and 
symmetric off-diagonal coupling 35 , and higher-order cou- 
plings. Work in this direction is now in progress. 



where the first two terms connected with the time deriva- 
tives can be calculated as follows 

(D 2 (t)\-\D 2 (t)) = 5>„i 2 ]Ta;a 9 



•El^| 2 E[-^(A,A;+A g A*)], 

n q 

(A2) 



and 



D 2 (t)\-\D 2 (t)) = E^-E^«l 2 E A '« A s 



Ei^i 2 E[4( a ^ + a ^)] 



(A3) 

The remaining term is the average energy in the D2 trial 
state: 

{D 2 \H\D 2 ) = (D 2 \H CX \D 2 ) + (D 2 \H ph \D 2 ) 

+ (D 2 \H cx _ ph \D 2 ) (A4) 



with 



(D 2 \H CX \D 2 ) = - J E VCVWl - ^E ^n^n-l 
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Appendix A: The D2 trial state 

The time evolution equations for Di, D and D2 
Ansatze can be derived by employing Dirac-Frenkel time- 
dependent variation method. We start from the simplest 
of the three, the D 2 Ansatz. For D2 Ansatz, the La- 
grangian is defined as 



L = (D 2 {t)\^-^-H\D 2 (t)) 

= \\(D2(t)\ ¥t \D 2 {t)) - (D 2 (t)\-\D 2 (t))] 
-(D 2 (t)\H\D 2 (t)) (Al) 



(D 2 \H ph \D 2 ) -E^«| 2 E W ?I A ?| 2 

n q 

{D 2 \H ex _ ph \D 2 ) = E l^| 2 E-9? w ?(V ,;<? " + K e ' iqn ) 

n q 

Equations of motion from the trial wave function D 2 (t) 
are readily obtained from Eq. (1231) 

+Jip n +l{t) + Jlp n -l(t) 

-1> n (*) E IK (*K 9 " + K (t)e~ iqn \ 

q 

(A5) 

-i\ q (t) = -ElV'nl'e-^^-^A, (A6) 

n 

It can easily be shown the norm of the D2 trial state is 
conserved, i.e., 

-§[f>n(*)| 2 ]=0 (A7) 
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And in this paper, we set 



JY 



(A8) 



Solutions to Eq. (|A5|) and (|A6j) provide the dynamic in- 
formation of the exciton-phonon system. 



Appendix B: The D trial state 

For the D trial state, the localized backbone of the 
translationally- invariant GL Ansatz, we can also derive 
coupled equation for the time-dependent variational pa- 
rameters using the Lagrangian formalism 



, ~ , , ,ih d 
L = (D(t)\-— - H\D(t)) 

= j{(m\^ t \D(t))-{D(t)\-\D(t))] 



-m)\H\D(t)) 



(Bl) 



wherein the individual terms can be calculated as follows 

- 1 - 
(D(t)\-\D(t)) 

n n q 

-a ? )(/3 9 v«"~a;)] + ]T|v>„| 2 

n 

9 

+(p q e-*^ - A)(/3*e^ - A,)] (B2) 



- V - 
(D(t)\-\D(t)) 



^iqn 



n n q 

-\*)(P q e-^-\ q )]+J2Wn\ 2 

n 

9 

+ {P q e- iqn - \W q e iqn - A,)] (B3) 



and 



(D\H\D) = (D\H CX \D) + {D\H ph \D) + (D\H cx - ph \D) 

(B4) 

with 

(D\H CX \D) = -J^C(i+iS n 



(D\H ph \D) = n- 1 e m 2 Y.^ e ' iqn - w 



(D\H e ^ ph \D) = iV- 1 /2^| V , n |2^ g?a;g[((3 * e i9n_ A * 



e' tqn + {/3 q e- iqn - A*)e 



Substituting Eqs. (|BTj) - ([B4|) into Eq. (|23)). one arrives 
at the equations of time evolution for the D Ansatz: 



-0*e^ - X*)(p q e-^ - X q )} 

+Jlpn+lS n ,n+l + Jlpn-lS n ,n—\ 

-N- 1 ^ n J2^\l3< ! e- iqn - X q \ 2 - N- 1 / 2 ^ 
i 

5> g w„[^ - K eign + V i9n ] ( B5 ) 



iN- 1 Y<\^\ 2 \{t) = 

n 

n 



—iqn _ p —iq(n+ 



+ -N- 1 jY,^n-lS n!n -l 



-iqn _ p -iq(n- 



1] ] 



+N- 1 Y,\^n\ 2 ^ q {P q e-' lqn - \) 

n 

+N- 1 / 2 Y / tt>n\ 2 g q u q e- iqn (B6) 

n 

n 

n 

+JN- 1 J2<i>n+lS; 



iqn 



n,n+l 



[^(e-^-l) + -A g e^(e^-l)] 

+ JN~ 1 J2^n-lS n , n -l 
n 

[P q {e lq - 1) + \\e^ n {e^ - 1)] 

-N- X Y^\^ q {Pq-Xqe iqn ) 

n 

-N- x l 2 Y,\i>n\ 2 9qUq (B7) 
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Appendix C: The Di trial state 



where 



The same procedure in Appendices E] and [B] can be 
applied to the Di Ansatz 36 with variational parameters 
ip n and X qn , and one can obtain 



+ Jlpn+lS n , n +l + Jlp n -iSn,n-l 

-N- 1 ' 2 ^ 9 q u q {\ qn e iqn + c.c. 



f(t n ,a(t n )) = 



da(t n ) 
dt n 



or A g , and is ob- 



and from t„ to t n +i = t n + At. 

For the D2 Ansatz, a represents ipn 
tained by Eqs. (|A5I) and (| A6|) . Then, by calculating ki, 
k 2 , k% and ki, one can get the parameters a(t n +\) for 
the next time step using the fourth-order Runge-Kutta 
method. 



(Cl) 



Appendix E: Linear absorption spectrum by D2 and 
D Ansatze 



-IN' 1 \lP n \ 2 \ qn {t) = 

+ N~ 1 Jlpnlpn+l(\q,n+l ~ \n)S n ,n+l 
+ A r_1 J-0*V'«-l(A g , n _l - X qn )S n ,71 — 1 

-N~ 1 \ip n \ 2 ui q \ qn 

-N-^l^g^e- 1 "" (C2) 

These coupled differential equations can be numeri- 
cally solved by the fourth-order Runge-Kutta method, 
and work on this is now in progress. 

Appendix D: Numerical Details 

There are various numerical approaches to solve cou- 
pled differential equations such as Eqs. (|A5j) and (|A6|) for 
the D 2 Ansatz and Eqs. ([B5|), (|B6} and JET} for the D 
Ansatz. One way is to transform these equations into the 
Volterra integral equations 3 ^ for minimization purposes, 
and then use a nonlinear optimization method such as 
the Newton-Raphson method to solve them. Another 
approach is to solve the time-dependent differential equa- 
tions directly. The latter method has a much higher com- 
putational efficiency, but it requires a higher precision in 
the single iterative time step since the computational er- 
ror may accumulate as the number of the iterative steps 
increases. 

In this paper, we use the Runge-Kutta fourth-order 
method 3 ^ to solve Eqs. (jMJ) and §M§ and Eqs. (|B5l) . 
(|B6|) and (|B7j) . The Runge-Kutta fourth-order method 
is widely used to solve differential equations. Its single 
step error is fifth order, i.e., 0(At 5 ). The algorithm for 
this method can be described as follows: 



ki 


= f(tn, 


a{t n ))At 






= f(t n - 


\--At,a(t n ) 


+ -*i)Ai 


h 


= f(t n - 


f -At,a(t n ) 


+ \k 2 )At 


fc 4 


= /(«»- 


\-At,a(t n ) + 


h)At 


a(t n+ i) 


= a{t n ) 


fci k 2 
+ — + — + 
6 3 


| + | + 0(A* 5 ) 



The linear absorption spectrum of the exciton-phonon 
system studied in this paper is calculated by 

F(u) = -Re / F(t)e lult dt (El) 

7T JO 

with 

F(t) = P h(0|ex(0|Pe-^ t P t |0) cx |0) ph (E2) 
where P is the polarization operator 

P = /i^(|7l)cx ex(0| + |0) O x ex(n|) (E3) 
n 

Here fj, is the transition dipole matrix element for a single 
site, and |n) cx is the exciton state at the nth site |n) e x = 

Substituting Eq. (|E3|) into Eq. (|E2l) one obtains 
F(t) =V 2 J2J2 ph(0|cx(m|e^ iAt |n)ex|0) ph (E4) 

n m 

wherein e _! - fft |n) cx |0) p h can be approximated by a Davy- 
dov trial state, for example, by a D2 trial state: 

e-^ t |n)cx|0) ph « 

y^n'-n(£)K)cx 
n' 

x exp(^[A 9 (t)6t - A* (*)£,]) |0} ph (E5) 
<? 

where the variational parameters ip n >- n {t) (n' — 
0, 2, N-l) and X q (t)(q = -N/2+1, -N/2+2, N/2) 
have the following initial values 

l(>n'-n{t = 0) = 5„/_„, (E6) 

and 

\ q (t = 0) = 0. (E7) 
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Substituting Eq. (IE5|) into Eq. (|E4|) . we obtain the for- 
mula to calculate F(t) by D 2 Ansatz: 

n m 

x ph (0|exp(^[A g (t)S+-A*(t)S ? ])|0) ph 
i 

n m q 

= ^Nj2^ m (t)exp(-^J2\ X M 2 ) (E8) 

m q 

wherein X q (t) and ip m (t) are initialized by Eq. (|E7j) 
and Eq. (|E6|) (i.e., Vm(0) = <^m,o) an d then solved by 
Eqs. and (jA"5l 

The same procedure from Eq. (|E5P to (|E8|) can be ap- 
plied to D Ansatz, and one obtains 

F(t) = >x 2 Nj2 iM*) ex P (-i E l^(*)^ i9m - A 9 (t)| 2 ) 

rn q 

(E9) 

where the variational parameters ?/>,„(£), /3 q (t) and A g (£) 
are solved by Eqs. ()B5|) - (|B7j) with the initial values 
iMO) = <W, /3 g (0) = and A g (0) - 0. 



Appendix F: Deviation vector of the D2 Ansatz 

Substituting Eqs. ([H)-(|3[l. (flS) and (HHJ) into Eq. ([52]) . 
one obtains the expression of |<5(£)) for the D2 Ansatz: 

\s(t)) = Y,B}MMt) + 

n 

+J[lp n+1 (t) +Vn-l(*)] - ?An(i)E W 9^ A ?W 

exp(^[A g (t)fet _ A;(i)6 g ])|0) ph |0) cx (Fl) 
9 
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